%this file is to calculate 2-point correlations.
%It does not consider homogeneity in any directino, all
%calculations are carried out via direct correlation.

%Created by Ricardo Mejia-Alvarez.
%University of Illinois at Urbana-Champaign
%Urbana, IL. 06/16/2010

clear

CL = 0.5;   %factor to find the correlation length

%input data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
display('how many different data folders are you going to consider?');
NF = input('?  ');
folData = cell(1,NF);
pathFiles = cell(1,NF);
NOR = zeros(1,NF);

for k = 1 : NF
   
    kstr = num2str(k);
    message = strcat('Enter folder number :',kstr,' \n');
    folData{k} = input( message , 's');
    
end


display('Now, for each data folder, tell me what columns to use to calculate')
display('two-point correlations (use vector format).')

for k = 1 : NF
   
    kstr = num2str(k);
    message = strcat('Columns for data folder number :',kstr,' \n');
    DataColumns{k} = input( message ); %#ok<*SAGROW>
    cols(k) = length(DataColumns{k}); 
end

cols = sum(cols);

disp('Enter the location of reference traces in cell format.')
Loc_RC = input('?');
tr = length(Loc_RC);

for tt = 1 : tr
    LRC_tr(tt) = length( Loc_RC{tt}(:,1) );
end



disp('Enter data transformation if necessary. Use matrix format.')
disp('For instance: ')
disp('                    X  Y  U1  U2  U3  ...  Un')
disp('Inversion/scale   [-1  1  -1   1   1  ...   1 ]')
disp('Horizontal flip   [ 1  0   1   0   0  ...   0 ]')
disp('Vertical flip     [ 0  1   0   0   0  ...   0 ]')
disp('Transpose         [ 0  1   0   0   0  ...   0 ]')
disp('  ')
disp('Enter only the cells to modify.')
disp('E.g. T(1,1)=-1; T(1,4)=-1; T(2,1)=1; T(2,4)=1; T(3,2)=1; T(4,2)=1')
disp('   ')
disp('Type "return" and press the Return key to finish')
T = ones(1 , cols + 2);
T(2 : 4 , :) = zeros(3 , cols + 2);

keyboard

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%creating the directory structure of the analysis
[~,folResults] = newfolder(folData{1},'results'); 


for k = 1 : NF

    %generating a vector with the name of the files to process
    pathFiles{k} = filesVector(folData{k},'*.dat');
    
    if numel(pathFiles{k}) == 1
        pathFiles{k} = filesVector(folData{k},'*.V3D');
    end

    %finding the number of instantaneous realizations
    NOR(k) = length(pathFiles{k});

end

testLength = sum(NOR - NOR(1));
if testLength ~= 0
    display('The number of files in the different folders is not the same.')
    display('I think you should take a look at this.')
    display('hit Return to continue')
    pause
    testLength = [];
end

NOR = NOR(1);



tic     %measures total time
tic     %measures calculation + data loading time


for tt = 1 : tr
    Location_RC = Loc_RC{tt};
    LRC = LRC_tr(tt);

    SigmaVel = [];
    for n = 1 : NOR

        p = 0;
        U = [];
        for k = 1 : NF
            
            filename = char(pathFiles{k}(n));
            [~,I,J,Dx,Dy,A] = matrix(filename);

            for m = 1 : length( DataColumns{k} )

                U{p + m} = A{ DataColumns{k}(m) };

            end
            p = length(U);

            if k == 1 && n == 1
                [A(1 : 2)] = Transformation(T(: , 1 : 2) , A(1 : 2));

                X = A{1};
                Y = A{2};
                
                DX = abs( X(1,2) - X(1,1) );
                DY = abs( Y(2,1) - Y(1,1) );

                X = repmat(X , [1 , 1 , LRC]);
                Y = repmat(Y , [1 , 1 , LRC]);

            end

        end
        
        [U] = Transformation(T(: , 3 : end) , U);

        for mm = 1 : LRC
            for nn = 1 : p
                V(mm , nn) = U{nn}( Location_RC(mm,1) , Location_RC(mm,2) );
            end
        end
       

        [Uc] = cross_variables(U , V);


        if n ~= 1

            for kk = 1 : LRC
                Ucorr{kk} = Ucorr{kk} + Uc{kk};
            end
            
            SigmaVel = SigmaVel + V.^2;

        else

            Ucorr = Uc;
            SigmaVel = V.^2;

        end

    end
    

    toc

    SigmaVel = sqrt( SigmaVel / NOR );
    Sigma1 = reshape(SigmaVel' , [p , 1 , length(Uc)]);
    Sigma1 = repmat( Sigma1 , [1 , p , 1] );
    Sigma2 = reshape(SigmaVel' , [1 , p , length(Uc)]);
    Sigma2 = repmat( Sigma2 , [p , 1 , 1] );

    SigmaVel = Sigma2 .* Sigma1;
    clear Sigma1 Sigma2
    SigmaVel = reshape( SigmaVel , [1 , p^2 , length(Uc)] );

    LcorrH = zeros(LRC , p^2 + 2);
    LcorrV = zeros(LRC , p^2 + 2);

    for kk = 1 : LRC

        Ucorr{kk} = Ucorr{kk} / NOR;
        
        LcorrH(kk , 1 : 2) = [ X(Location_RC(kk , 1) , Location_RC(kk , 2)) ,...
                Y(Location_RC(kk , 1) , Location_RC(kk , 2))];
            
        LcorrV(kk , 1 : 2) = [ X(Location_RC(kk , 1) , Location_RC(kk , 2)) ,...
                Y(Location_RC(kk , 1) , Location_RC(kk , 2))];

        for jj = 1 : p^2
            Ucorr{kk}(: , : , jj) = Ucorr{kk}(: , : , jj) / SigmaVel(: , jj , kk);
            
            CL1 = CL * Ucorr{kk}( Location_RC(kk , 1) , Location_RC(kk , 2) , jj);
            
            ML = Ucorr{kk}( Location_RC(kk , 1) , 1 : Location_RC(kk , 2) , jj);
            MR = Ucorr{kk}( Location_RC(kk , 1) , Location_RC(kk , 2) : end , jj);
            MU = Ucorr{kk}( 1 : Location_RC(kk , 1) , Location_RC(kk , 2) , jj);
            MD = Ucorr{kk}( Location_RC(kk , 1) : end , Location_RC(kk , 2) , jj);

            LL = find( abs(ML - CL1) - min( abs(ML - CL1) ) == 0 , 1 , 'last' );
            LR = find( abs(MR - CL1) - min( abs(MR - CL1) ) == 0 , 1 , 'first' );
            LU = find( abs(MU - CL1) - min( abs(MU - CL1) ) == 0 , 1 , 'last' );
            LD = find( abs(MD - CL1) - min( abs(MD - CL1) ) == 0 , 1 , 'first' );
            
            LcorrH(kk , 2 + jj) = (length(ML) - LL + LR - 1) * DX;
            LcorrV(kk , 2 + jj) = (length(MU) - LU + LD - 1) * DY;
        end

    end



    toc

    tic     %measures saving-time

    disp('saving results')

    ti = num2str(I);
    tj = num2str(J);

    LocStartX = num2str( Location_RC(1 , 1) );
    LocStartY = num2str( Location_RC(1 , 2) );
    LocEndX = num2str( Location_RC(end , 1) );
    LocEndY = num2str( Location_RC(end , 2) );


    filename = strcat('\Non-Homog_[',LocStartX,',',LocStartY,']_to_[',LocEndX,',',LocEndY,'].dat');

    for k = 1 : LRC
        zoneX = num2str(Location_RC(k , 1));
        zoneY = num2str(Location_RC(k , 2));

        heading = strcat('ZONE T = "[',zoneX,',' , zoneY ,']" I = ' , ti , ' J = ', tj , ' F = POINT');

        data = cat( 3 , X(: , : , k) , Y(: , : , k) , Ucorr{k} );

        data = mixingMat(I , J , data);
        data = sortrows(data , [2 , 1]);

        if k ~= 1
            saverappend(folResults,filename,heading,data)
        else
            saver(folResults,filename,heading,data)
        end

    end
    
    ti = num2str(LRC);
    heading = strcat('ZONE I = ' , ti , ' F = POINT');
    
    filename = strcat('\CorrLengthHorz_[',LocStartX,',',LocStartY,']_to_[',LocEndX,',',LocEndY,'].dat');
    saver(folResults,filename,heading,LcorrH)
    
    filename = strcat('\CorrLengthVert_[',LocStartX,',',LocStartY,']_to_[',LocEndX,',',LocEndY,'].dat');
    saver(folResults,filename,heading,LcorrV)
    
    clear LcorrV LcorrH Ucorr Uc V X Y A data
    
    save
    clear all
    load

end
    toc
    
